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ABSTRACT 


In this thesis we derive a lattice structure to realize linear phase transfer 
functions and develop an adaptive algorithm for continuously updating the lattice 
reflection coefficients. The lattice structure is considered because of its superior finite 
wordlength performance compared to transversal structures. The adaptive lattice 
algorithm developed in this thesis has been applied to estimate the sinusoidal 
frequencies as part of Prony’s method. Results of computer simulation supporting the 


theory are reported. 
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Il. INTRODUCTION 


Every finite-impulse response (FIR) filter has two distinct properties [Ref. 3] ; 
first, it is always stable; second, if it is not causal, it can always be made to be causal 
by introducing a finite delay. FIR filters can be designed so that their frequency 
responses have an exactly linear phase characteristic. FIR filters with linear phase are 
important in applications like speech processing and data transmission, because in 
these applications a nonlinear phase filter is harmful. The symmetry property of the 
linear phase FIR filters, however, helps reduce the number of coefficients by nearly one 
half resulting in substantial computational savings. 

The various filter realizations, or structures that are frequently considered are the 
direct form, the cascade form, the parallel form, and the lattice form. The lattice form 
realization is of particular interest because of its superior numerical performance, and 
modularity in the structure. The operation of a lattice filter is completely described by 
specifving the sequence of reflection coefficients that characterize the individual stages 
of the filter. 

Conventionally an adaptive filter is composed of a tapped delay line or 
transversal structure with adjustable coefficients or weights and an adaptive algorithm 
which updates the coefficients continuously based on some performance criterion. The 
design of a fixed coefficient filter is based on the prior knowledge of both signal and 
noise. Adaptive filters, on the other hand, have the ability to adjust their own 
parameters automatically, and their design requires little, or no a priori knowledge of 
signal or noise characteristics [Ref 14]. However, the designer has to choose the order 
of the filter and the type of the algorithm. Also the adaptive filter usually requires a 
large initial transient time (1.e., the initial filter convergence period). 

The least-mean-square (LMS) adaptive algorithm minimizes the mean square 


ry 


error £(k) by recursively altering the filter coefiicient vector 3(K}) at each sample 
instant. The original Widrow-rioff LMS algorithm is B(k+ 1)=8(K) + 2pek)A(k), 
where [t is a convergence factor controlling stability and rate of adaptation [Ref. 15]. 
The algorithm is based on the method of steepest descent, moving 6(k) in proportion 
to the instantaneous gradient estimate of the mean square error. The successive 
orthogonalization provided by the lattice offers advantages in adaptive convergence 


rate which cannot be achieved with tapped delay lines [Ref. 17,18]. 
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Recently, Prony’s method has been applied to the estimation of spectral lines in 
noise [Ref. 11]. It has been shown that the Prony’s method yields better spectral 
estimates than a companion spectral estimation technique, called Pisarenko’s method 
[Ref. 10]. Also it has been established that the filter structure involved in Prony’s 
method has a linear phase property which is not the case for Pisarenxo’s method [Ref. 
11}. 

The thesis investigates the application of lattice structures in Prony’s method of 
spectral line estimation. The complete solution to Prony’s method consists of three 
steps: (1) representing a given process of M sinusoids in noise in terms of complex 
exponentials, (11) finding roots of a symmetric polynomial, and (i) estimating the 
frequency, phase and amplitude information. In this thesis, however, we have 
emphasized only the frequency estimation problem. Here we derive an adaptive lattice 
structure to realize linear phase transfer functions which will be used to estimate the 
sinusoidal frequencies as part of Prony’s method. The scope of the thesis consists of 
obtaining a linear phase lattice structure, developing a least mean square (LMS) error 
based adaptive algorithm, and testing the lattice structure and the algorithm by means 
of computer simulation. 

The thesis is organized as follows. In Chapter II, we present a brief review of 
Prony’s method of representing a given process in terms of a set of complex 
exponentials, and then address the problem of estimating spectral lines using this 
method. We show that the original Prony’s method has to be modified slightly in 
order to apply it to the spectral line estimation problems. In Chapter III, we discuss 
the basic concepts of the linear phase FIR filter and the related lattice structure. We 
show three examples of obtaining the lattice structures for both linear phase and non- 
linear phase FIR transfer functions (Appendix A). In Chapter IV, we briefly discuss 
the least-mean-square (LMS) adaptive algorithm that results from attempting to 
minimize the mean square error and present a summary of some LMS algorithms for 
lattice that have been reported in the literature. Also included are some computer 
simulation results. [n Chapter V, we present the derivation of a new LMS based FIR 
adaptive lattice aigorithm and extend it to the linear phase case as weil. ihe new 
algorithm is then used in the estimation of spectral lines in white noise. Results of 


simulation are included. 


i 


Ii. PRONY’S SPECTRAL LINE ESTIMATION 


A. INTRODUCTION 

In its original form, Prony’s method analyzes processes involving simple, or 
damped harmonics using complex exponential functions as coordinate functions. On 
the other hand, the well known Fourier analysis consists of representing a given 
process in terms of a set of sine and cosine functions. Recently, Prony’s method has 
been applied to the estimation of spectral lines in noise [Ref. 11]. [It has been shown 
that the Prony’s method yields better spectral estimates than a companion spectral 
estimation technique, called Pisarenko’s method, in terms of bias, spurious responses, 
and the computational complexity [{Ref. 10]. 

In this chapter, we present a brief review of the Prony’s method of representing a 
given process in terms of a set of complex exponentials, and then address the problem 
of estimating spectral lines using this method. We show that original Prony’s method 


has to be modified slightly in order to apply it to the spectral line estimation problems. 


B. PRONY’S METHOD 
Consider that the given process, x(k), consists of n distinct sinusoids, then x(k) 


can be approximated by an expression of the form 


n 
x(k) =, 2% [A, cos@.k + B. sinw.k] (25 
i 


pores 


where the @,’s are sinusoidal frequencies. The above approximation can be considered 


as a special case of an exponential approximation given by [Ref. 29] 


im a Ne 
x(k) = SC. els (2.2) 


i=] ! 


—" 


where m = 2n and the a.’S are identified as oO, and -O.. The values of a, can be 
estimated, by Prony’s method, assuming that the data are known at k = 0, 1,..., N-l. 


From Eq.(2.2), we can obtain the following set of equalities 
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C, el) + C, e242 + .. + om ean = ich) 


meee] + Ci cl*49 +. +C_ el@4m = = (2) (2.3) 


C, elN-Day +c, el(N-Day 4 +c eMN Dan = (N-1) 


Now we have a set of N equations with 2m unknowns, namely, CeenGe a Game ns 
m) which can be solved exactly if N=2m, or approximately by the method of least 
squares if N>2m. Also note that the N equations are nonlinear in the exponential 


terms e!4i, Let el, lowly 2.4 ta, OC tae roots Of the equation [Ref. 29}. 


elma fe g,eKm-l)a ate ee Stee Sic caie a y.je!" a Gon = 0 (2.4) 


In order to determine the coefficients a. (i = 1, 2,..., m), we multiply the first equation 


ieteG.(2.3) by a» the second equation by eetne mth equation by a, and the 


m-1? 
(m+ yyth equation by 1, and add the results. In this way, we can obtain the (N-m) 
linear equations. From Eq.{2.2), we can obtain the following set of equalities 

x{m) + x(m-1)a@, + x{(m-2)a, + .. +x(O)a, = 0 


x(m+i) + x(m)a, + x(m-I)a, + .. +x(I)a, = 0 


(2.5) 


x(N-1) + x(N-2) a, + x(N-3) a, + .. +x(N-m-l)a@, = 0 


Since the data x(k), k = 0, 1, 2,..., N-l, are known, this set of equations can be solved 
for them a@’sif N22m. After the a’s are determined, el, 1 = I, 2,..., m, are found 
from Eq.(2.4). The set of equations in Eq.(2.3) then becomes linear in terms of C.,i = 
1, 2, .., m. The C’s can be determined from the first m ecauations, Oredeteriusee 
approximately by appiying the least squares method to the entire set of equations. The 
nonlinearity associated in finding el, which is related to the frequencies j@, and -jq,, is 
concentrated in Eq.(2.4). The above procedure is known as Prony’s method. 

Now, in the case of Eq.(2.1), since the roots of Eq.(2.4) occur in reciprocal pairs, 
Eq.(2.4) must be invariant under the substitution of eJ4: for el3;, so that We must have 


Uae Hn) = Ops er Ong 1 = Eey: Thus Eq.(2.4) becomes 


el2no + chy el(2n-O 1 4 Gn] aint Io + an eine 
ene en-HO 4 4 ay J) +1 = 0 


OT 


QIN f (eINM 4 IMM 4 g, (ellO DO + ela) + 4+ 
Gy (el eet ae 


since J2® + 0 we have 


2cosn® + 2@, cos(n-ljw +... + 2 | COS@an Cameo (2.6) 
Now noting that m= 2n, and applying the above symmetry of a@’s to Eq.(2.5) yields 
UNG 2 cae 
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{x(0) + x(2n)} + {x(1) + x(2n-])} a, +... + {x(n-l) + x(n+])} ay 


mend, = 0 


{x(1) + x(2n+1)} + {x(2) + x(2n)} @, +... + {x(n) + x(n+2)} a,_, 


ment ij)e, = 0 


eo 


{x(N-2n-1) + x(N-1)} + (x(N-2n) + x(N-2)} a) +. + {x(N-n-2) 


meen); @. tf x(N-n-l)a, = 0 


Eq.(2.7) consists of a set of N-2n equations and n unknowns. This set can be solved 
exactly for the a’s if N=3n, or solved approximately by the least squares method if 
N>3n, and then the @’s are determined from Eq.(2.6). 


C. ESTIMATION OF FREQUENCY, AMPLITUDE AND PHASE 

If a process under measurement contains an unknown number of sinusoids of 
unknown frequencies and amplitudes, a variant of Prony’s method can be used to 
determine the number of sinusoids and their associated frequencies and amplitudes. As 
noted above, Prony’s method is a technique for modeling data of equally spaced 
samples by a linear combination of exponentials. An exponential curve having M 
exponentials terms can be determined from the 2M data measurements. Each 
exponential term A.eAK has two parameters —~ an amplitude A, and an exponent a, 
where a, can be real or imaginary. For the case where only an approximate fit with M 
exponentials to a data set of N samples is desired, such that N>2M, a least squares 


estimation procedure is used. 


The model assumed is a set of M exponentials of arbitrary amplitude, phase 


frequency and damping factor. A process consisting of M real undamped (a is 
imaginary) sinusoids can be expressed as 


x(k) = . ¥ (az +a, 7,5) (2.8) 


—_ A, cos(2m£kT + 0.) 


with @. = (A,/2)e!% and Z,= Ql2rt T where A, is the amplitude, f, is the frequency and 0 
is the aa of the ith sinusoid, respectively, and T is the sampling interval 

Finding {A., 0 f.} and M that minimize the squared error is a difficult nonlinear 
least squares problem. Prony’s method solves two sequential sets of linear equations 


with an intermediate polynomial rooting step that concentrates the nonlinearity of the 
problem in the polynomial rooting procedure (similar to Eq.(2.4)) 


Define the polynomual, B(z), which has z, and Ze as its roots, given by 


B(z) = i (2- 2)(2- 2, “) = ¥ 5° 72M-j = 9 (2.9) 


with bo = 1 and the b, being real coefficients. Since the roots are of unit modulus occur 


in complex estan pairs, then Eq.(2.9) must be invariant under the substitution z™ 


-{ 
for z. Therefore, Eq.(2.9) can be written as 


2M 2M 
22M Bi /z) = 22M £ bz 72M = a” b, 2! = 0 (2.10) 
ene a 


Comparing Eqs.(2.9) and (2.10), we may conclude that b= bou-j for }=0,1, <2 aye 
with by = bo,, = |. Thus the requirement for complex conjugate root pairs of unit 


modulus is implemented by constraining the polynomial coefficients to be symmetric 
about the center element. 


Hence the realization of B(z) is a linear phase FIR filter 
Based on order 2M, a linear prediction error can be written as 
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V 
e(k) = Eb, [x(k +j) + x(2M-K+)) 


Cet 
— 


ay 


which reduces the number of coefficients required by one-half. 


The coefficients b,, b,, .. . 6 are determined in a least squares fashion by 


minimizing the total squared error. 


N-2M_ 5 7 
E= 2) eth) (2.12) 
|= 


which yields the normal equations 


N-2M 


, yl ke [X2M-k+i) + x(k +i] [X2Mj+i) + xj+H] = 0 (2.13) 
E el 


for k=1,... ,M 


This equation can be solved recursively for increasing order M. 


From the estimated {b;} values, the {z,} are determined using Eq.(2.9). This gives 
the frequency estimates. 


7 tan} [Im(z,)/Re(z,)] / 2nT 


(2.14) 
To obtain the {d. } a second set of normal equations is solved, 
; x [ a. x z,1 zJ + a = za a 
i=l 'j=0 * |} /  7=0 
N CO 
= 2 (2 Rez,J) x(j) (ats) 
ia” 


fOr k= 0.105... Vi 


The set {0.} provides both amplitude (A,),or power, and phase (0.) information: 


fa el |e | (2.16) 

@. = tan"! [ Im(b,)/Re(b,) | (2.17) 
In the foregoing we have shown that the frequency estimation as part of Prony’s 
method requires a polynomial with a linear phase property. In the next chapter, we 


show that a lattice structure can be utilized to implement a linear phase transfer 


function. 
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lil. LINEAR PHASE FIR FILTERING 


A. INTRODUCTION 

Every finite-impulse response (FIR) filter has two distinct properties {Ref. 5] ; 
first, it 1s always stable; second, if it is not causal, it can always be made to be causal 
by introducing a finite delay. 

FIR filters can be designed so that their frequency responses have an exactly 
linear phase characteristic. FIR filters with linear phase are important in applications 
like speech processing and data transmission, because in these applications a nonlinear 
phase filter is harmful. 

The impulse response sequence of a linear phase FIR filter exhibits a kind of 
symmetry, e.g., h(n) = h(N-1I-n) for n = 0, I, ..., (N/2)-1 (assuming that N is even). 
In general, FIR filters require a large number of coefficients to adequately meet with 
tae required filter specifications. The symmetry property of the linear phase FIR 
filters, however, helps reduce the number of coefficients by nearly one half resulting in 
substantial computational savings. 

The various filter realizations, or structures that are frequently considered are the 
direct form, the cascade form, the parallel form, and the lattice form. The lattice form 
realization is of particular interest because of its superior numerical performance, and 
modularity in the structure. Lattice realizations have been successfully utilized in 
filtering and spectral analysis, and in modeling of some physical process like speech, 
geophysical data etc. [Ref. 8]. The operation of a lattice filter is completely described 
by specifying the sequence of reflection coefficients that characterize the individual 
stages of the filter. 

In this chapter, we will discuss the basic concepts of the linear phase FIR filter 
and the related lattice structure. And finally, we will show three examples of obtaining 
the lattice structures realizing for both phase and non-linear phase FIR transier 


functions (Appendix A). 


B. LINEAR PHASE FIR FILTERS [REF. 3] 
Let {h(n)} be a causal finite duration sequence defined over the interval 
OSn<N-I1I having the linear phase symmetry property given by h(n) = h(N-I-n). The 


z transform of {h(n)} can be written as 


le 


N-1 
2 


H(z) = d(njz™ end) 


if N is even, then Eq.(3.1) can be rewritten as 


(N/2)-1 : N-1 ; 
EZ) =e) ee nye DY Ah un 
(Z) ie (n)z = (n)z 


Now applying the symmetry property of h(n) in the second term on the right hand side 
yields 


N/2)-1 N/2)-1 
H(z) Ig h(n)z2 NY) h(N-I-n)z(N--n) 
n=0 n=0 
which can be simplified to 
(Ni2)- 


H(z) = = - h(n) {22 + z(N-1-n) (3.2) 


If N is odd, then Eq.(3.1) becomes 


H(z) = UN-/2h | h(n)[iz + z(N-}-n)y + hi) zI(WN-D/2] (3.3) 
es 


If we evaluate Eqs.(3.2) and (3.3) for z= J , we obtain the discrete Fourier transform 


of the filter sequence h(n) defined as 


Sei 
H(el®) = E h(n) eJOn (3.4) 
n= 


For the case when N 1s even, we then have the discrete Fourier transform 
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Baie Lyi \ N/?)- 
H(el@)=e ee _. 2h(n)cos{ o{n-(N-1)/2)}]] i ae)) 


and for N odd, we obtain 


vay) = a JM{(N-1)/2 N-3)/2 
H(el®) = &7 UND ENT V2} 4 | e 2h(n) cos{oo{n-(N-1)/2}]}. (3.6) 
n= 

In both cases above, the sums in brackets are real, implying a linear phase shift 
corresponding to a delay of (N-1)/2 samples. Figures 3.1 and 3.2 show the direct form 


realizations of an FIR filter with linear phase for both N even and N odd cases. 
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Figure 3.1 Djirect-form realization for an FIR filter with linear phase (N= even). 
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Figure 3.2. Direct-form realization for an FIR filter with linear phase (N = odd). 


C. LATTICE FILTERS [REF. 4] 
The basic single section lattice structure is shown in Figure 3.3 where x(n) ts the 
x(n) and f,(n) and g)(n) are generally known as the forward and backward prediction 


errors, respectively. The defining equations of the lattice are given by 


f(n) = g(n) = x(n) 
f,(n) = fy(n) + Ky go(n-1) Ga) 


as kK fg(m) + g,(n-1) 


where Ky is called the lattice reflection coefficient. If another section is cascaded with 


the first one the resulting equations for the next order prediction errors are given by 
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pon) = f,(n) + kg, (n-1) 
g.(n) =e kof, (n) 1 g, (n-1) 


Substituting for f,(n) and g,(n-1) from Eq.(3.7) yields 
f(n) = x(n) + {k, +k, k,}x(n-1) + k,x(n-2) (3.8) 


or 
Pan) = x(n) + by 1X(n-1) ae bo 9x(n-2) 


and 
(eo) 


g.(n)=k,x(n)+ {k, +k,k,}x(n-1)+ x(n-2) 


or 
g.(n) = by 2x(n) a bo yx(n-1) + x(n-2) 


where b5 = Ky(1 +k») and b> q =k». 


Thus, by induction, for a cascade connection of M lattice sections we have the general 


expressions 


L by ;X(n-1) 
(G20) 


M 

— OMe) 
where by, p= 1. Eq.(3.10) represents the FIR filter type equations to obtain the yth 
Now taking the z-transform of 


order forward and backward prediction errors. 


Eq.(3.10) yields 
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MI : 
Pyiz) = * by iz  X(z) (330) 


{ : 
Gyy(Z) = = ONTn (3.12) 


=I 


For an unit unpulse input, Le., when A(z)= 1, Fy4(z) and Gy4(z) represent the forward 


and backward prediction error filter transfer functions, respectively, and 


Gy(z) = 2M Fah (3.13) 
eae) fem 
Fm (2) = Fm+ 1 (2) 
K+] 
Doe 
Gq (Z) Db v reine 
ra gm(k-1) Gi 





Figure 3.3" Uatticassection 


By inspection of Figure 3.3 we see that the Mth stage lattice reflection coefficient 
is equal to the FIR filter coefficient by, yg in Eqs.(3.11) and (3.12). In other words, 


we have 


24 


In fact, the FIR like prediction error filter coefficients can be iteratively calculated 
starting from Eq.(3.14). The algorithm, to calculate the filter coefficients BV | v also 
makes use of the weil known lattice property that an Mth order lattice contains all 
prediction filters of orders mS M. Now consider the m'} section (1SmsM) in the 


cascade connection of M lattice sections which can be described by the following 


equations: 
Fin(Z) = Fon.3(2) + ke !Grn_1 (2) (3.15) 
Gin(2) = KenFmn.p(2) + 27 Gey. 1(2) (3.16) 


Eq.(3.16) can be rewritten as 

era) — 20,{2) - zk,F | (2) Ca 
Substituting Eq.(3.17) into Eq.(3.15) yields 

eee 2) KG) kK a F -1 (2)} (3.18) 


Therefore, the (m-1)th order forward prediction error transfer function can be written 


th 


in terms of the m™ order forward and backward prediction error transfer functions as 


follows: 
eee) = KG (Zz) | 
F n-1 (2) se SS. (3.19) 
m 
where ea rem Pe ecaiiing 2 q.(3-13), that is, 
=i - | 
eZ Zz F(z *) (3.20) 


By substituting Eq.(3.20) into Eq.(3.19) we find that 
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F(Z) = Kazaa zh ) 


E (Z) = 
m-] 2 
ae Kem 


(3.21) 


Thus, from Eq.(3.21) we see that the next lower order polynomial F(Z) Cane 
calculated xXnowing F,,(z). Following the foregoing procedure we can find 
kena = Din- 1 m- 1; from F,,,_)(Z), and also obtain F,,,_5(z). This way, for a given yth 
order polynomial F,4(z), we can determine all the lattice reflection coefficients k_,, 
m=M, M-l, ..., 2, 1. This procedure is known as the step-down procedure [Ref. 1], 
and can be summarized as follows: Let the given M th order polynomial be 


th 


and by replacing M with m we have an m*™ order the polynomial for m lattice sections 


F(Z) = z Bae (3.23) 


where m= M, M-l,..., 2, 1. We now define another polynomial given by 


F(z!) = . bmn iz! (3.24) 


AS we Step through the procedure from m= M to m=1 Eq.(3.24) can be expressed as 


—— mae : 2 
ae oi = bnmiz ors 
Define k= bm and obtain the m-1' order polynomial as follows: 
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P(g) 2 ™e_ (z)) 


FF 4(z) = (326) 
m-l 2 
ek, 
Substituting Eqs.(3.23) and (3.24) into Eq.(3.26) yields 
“1 oo. 5-M m-i 
Pte ge Son oni 
3a ao 
|= } > 
ek, 
-i “j 
] Pm i2 ; kn; 2 Pm? 
“1 _ 
MF Pm = (3.27) 
ee hs 


m 


Then, by equating the coefficients of like powers of z on both sides of Eq.(3.27), we 


pth 


obtain the computational expression for the (m- order polynomual coefficients as 


b 





Det Da an 
el i wa a (3.28) 
Wo kA 


where m = M, M-l,..., 1 andi = 0, 1, ..., m-l, bo = enter and elec for a 
minimum phase polynomial F(Z). 


D. LINEAR PHASE FIR LATTICE FILTER 

In the foregoing we considered obtaining a lattice structure from a given FIR 
transfer function, and vice versa. In what follows we will deal with a special case of 
FIR filtering, namely, the linear phase FIR filter. Let {h({n)} be a causal finite duration 
linear phase sequence defined over the interval OSn<SN-l. From Eqs.(3.2) and (3.3), 


the z transform of {h(n)} can be written as 


Hy.(2) = Fy(z) + 2? Gylz) (3.29) 


cae 





Figure 3.4 Linear Phase FIR Lattice Fliter. 


Where Fa4(z) and G),(z) are forward and backward filter transfer functions, 
respectively, and D represents the number of unit delays. We now consider the N odd 
and N even cases separately. 

N odd case: For N odd, Eq.(3.29) can be written as 


Hy (2) = ag + ayz! + agz? +... + ani.3yg2 3? e ane pyz ee? 
_ angyaZ 1/2 4+ 0 +t ee) fe agai L agz (ND 


= ay + ay zi + az + te + ajy.3y/2 NP? + 2 (2) aN.pyge yD? 


ie acn.3)/22 ON * ee + az \N-3) a aie sta ak ee 
/ 
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Note that we are splitting the coefficient a(N-1)/2 into two coefficients of value 
(1/2)arn.1)/2 each. 


H~:_1{Z) = an a ai oT are mie tes ot arn.3y2z NP? foe) 


-+- 


(1/2) ay. pyg2 DP? + 2D! (172). 1/2 


+ 


aN.3/2Z tee aga + Bea 


+- 


ya) 1)/2 | 


Hy (2) = Fyg(z) + oN? Gy (2) (3.31) 


where the order of the polynomials Fy (Z) and Gyy(z), M = (N-1)/2, and the 
number of unit delays, D = (N-1)/2, the forward transfer function, 


2) = ag ae az * pce cso rhe an.3ygz ee? at (1/2)aqy. pyg2z 0D! (3:32) 
F(z!) a ay at az eral acn.3y/220 9)? a (1/2)arn. yD? (3433) 


and the backward transfer function, 


ss — + -(N-i)}/2 -1y Seah eee a 47k fey 
Gyj{Z) a ad Fy tz = Lan 1/2 BiNeoioe (5.04) 


sao Ie AAC) one 


Zo 


N even case: For N even, Eq.(3.29) can be written as 


= = s(IN- SING 
Ne) ae ag ae ayZ | ae agz - Sees rete SONED) 2a (N 2)/2 ae “(Ne ee N/2 


bn # agetN-3) + a2) 2 age NN 


: z af No2\/ 
Fins (2) = ag ae ayZ | alg az - oe baie SONED 2 (N 2)/2 (3735) 
ZENIT =f Ne ” 
ate N/2 [8(N-2)/2 eee e ots ayZ (\ 6)/2 
ea N-DIZ & age(N-2)2| 
Hy 1(2) = Fyyj(z) + 2° * Gy.(2) (3.36) 


where the order of the polynomials F,;(z) and Gyj4(z), M = ((N/2)-1}, and the 
number of unit delays, D = N/2, the forward transfer function, 


(N-2)/2 (3.37) 


= ° -2 
Fy 4 (Z) a ay = ayZ ate az ae aeons SONE2) 2a 


=) 
Fy 4(z Ly = ay a ayZ ar onze Tee. og a(n.2y/22™ 2)/2 (3.38) 


and the backward transfer function. 


Gy) = 2ON? Brel) = anayg + + agziN OW? (3.39) 


=P Re ae aie tagz(N-2)/2 
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Figure 3.4 shows a linear phase FIR lattice realization where the filter output is 
obtained by adding the Mth order forward prediction error, and the M‘® order delayed 
(oy D unit delays) backward prediction error. Combining the discussion in this section 
and that in Section (III.C), we can summarize the algorithm for converting a given 
FIR transfer function into a lattice structure as follows: 

femerind if N 1s even, or odd 

(i) For N odd: M = (N-1)/2, D = (N-1)/2 

(iii) For N even: M = {(N/2)-1}, D = N/2 

(iv) From Fy,4(z), obtain M reflection coefficients as discussed in Eqs.(3.22) to (3.26) 


(v) Implement the lattice as shown in Figure 3.4 


E. LATTICE REALIZATION OF A GENERAL FIR TRANSFER FUNCTION 

In the foregoing we have considered a polynomial F),(z) with buo= . However, in 
general we have bMo* 1. In this section, we shall modify the lattice realization 
algorithm presented in Section (III.C) to suit an arbitrary FIR transfer function with 
buo7 ! fcr 22|) The mA order polynomials F(z) and G,,(z) can be obtained 
mom 2 s.(3.13) and (3.16): 


F(Z) = sFy)(z) + Kin?” Gn. 1(2) (3.40) 
eee Km-1\2) + Sm 2 Om: {2) (3.41) 
where the coefficients s_. = bm,0 agyal Ns es Om,m are recognized as the reflection 

coefficients. Eq.(3.41) can be rewritten as 

(2) = Se 72 Ge 4) ie 5a ber ate) (3.42) 
Substituting Eq.(3.42) into Eq.(3.40) yields 

meee (Zz) +S” Kf Galz) - KF -1(2)} (3.43) 


Memererore, the (m-i)" order forward prediction error transfer function can oe 


th 


written in terms of the m™ order forward and backward prediction error transfer 


functions as follows: 


a1 





Figure 3.5 FIR Lattice. 


Sm En(2) = Kn Omi?) (3.44) 


Fin-1(2) = ae 2 
see k 


where s,,#k,, Substituting Eqs.(3.20), (3.23) and (3.25) into Eq.(3.44) yields 


- -m m-1 
: Sa Pm?" - kin? ya piyitiete 

“1 
"Soni? _ 


n ; 1 
=] -] 
“J _ om Y Pm, - km gee 
“1 
MF Ome FS (3.45) 
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From Eq.(3.45), we obtain the computational expression for the (m-1)th order 
polynomial coefficients is 
Semone = KD ; 
emai mm mt 2 16) 
Om-1.i a 5? : k? (5.46) 
m. m 
where m = M, M-l,..., 1 andi = 0, 1,..., m-l, os ores penne Onn sen Ae. 
for a minimum phase polynomial F,,(z). It may be noted that s,, = i form = I, 2, 


.., M-1. This indicates that we have only one s-coefficient, i.e., S\,, which requires an 


extra multiplication in the M*® Jattice section as shown in Pere sno: 


Aye 


IV. LMS ALGORITHM 


A. INTRODUCTION 

Conventionally an adaptive filter is composed of a tapped delay line or 
transversal structure with adjustable coefficients or weights and an adaptive algorithm 
which updates the coefficients continuously based on some performance criterion. 

The design of a fixed coefficient filter is based on the prior knowledge of both 
signal and noise. Adaptive filters, on the other hand, have the ability to adjust their 
own parameters automatically, and their design requires little, or no a priori knowledge 
of signal or noise characteristics (Ref: 14]. However, the designer has to choose the 
order of the filter and the type of the algorithm. Also the adaptive filter usually 
requires a large initial transient time (i.e., the initial filter convergence period). 

For stationary stochastic inputs, the mean square error, the difference between 
the filter output and an externally supplied input called the “desired response”. is a 
quadratic function of the filter coefficients, a paraboloid with a single fixed minimum 
point that can be sought by gradient techniques [Ref. 16]. 

In the previous chapter we showed that the operation of a multistage lattice filter 
is completely described by specifving the sequence of reflection coefficients that 
characterize the individual stages of the filter. In this chapter we briefly discuss the 
least-mean-square (LMS) adaptive algorithm that results from attempting to minimize 
the mean square error and present a summary of some LMS algorithms for lattice that 
have been reported in the literature. Also included are the computer simulation results. 

The least-mean-square (LMS) adaptive algorithm minimizes the mean square 
error €(k) by recursively altering the filter coefficient vector B(k) at each sampling 
instant. The original Widrow-Hoff LMS algorithm is B(kK+1)=8(k) + 2pte(k)X(k), 
where pH is a convergence factor controlling stability and rate of adaptation [Ref. 15]. 
The algorithm is based on the method of steepest descent. moving 6(k) in proportion 
to the instantaneous gradient estimate of the mean square error. 

When the filter input is stationary, the backward prediction errors are orthogonal 
to each other, with the result that the successive stages of the lattice filter are 
decoupled from each other [Ref 8]. This means that the global optimuzation of a 


multistage lattice filter may indeed be accomplished as a sequence of local optimization 


problems, one at each stage of the lattice filter. Accordingly, it is a straightforward 
matter to increase the order of the lattice filter by simply adding one or more stages 
Without affecting the earlier design computations. The successive orthogonalization 
provided by the lattice offers advantages in adaptive convergence rate which cannot be 
achieved with tapped delay lines [Ref. 17,18]. 


B. SUMMARY OF THE LMS ALGORITHM 

The LMS algorithm uses an estimate of the gradient of the mean square error obtained 
from the adaptive linear combiner which is a combination of a transversal structure 
and an adder. The adaptive linear combiner can be shown in two basic ways, 
depending on whether the input is available in parallel form (multiple inputs), or in 
serial form (single input) [Ref. 6]. 


d (k) 





4. y (k) e (k) 





Figure 4.1 Parallel Form. 


BD 





Figure 4.2 Serial Form. 


In the following we present a brief derivation of the LMS algorithm. Let us choose the 


single input form, then the filter output is given by 


y(k) = bo(k)x(k) + by(K)x(K-1) +. bye y(K)x(K-N+ 1) 


N-I | 
Ee by(k) x(K-i) 


X!(k) B(k) 


B'(KX(K) (4.1) 
where X and B are the input signal vector and the filter coefficient vector, respectively 


and (N-1) is the order of the filter. The input signal vector X and the filter coefficient 


vector B are defined as 
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x9) 1) Bok 
B(k) = 


X(k) = 


X(N +1) | by 1(k) 


Therefore, the output y(k) is equal to the inner product of X(k) and B(k). 
The error e(k) is defined as the difference between the desired response d(k) and 


the actual response v(x), 


e(k) = d(k) - X#(K)B(k) = d(k) - BE(K)X(k) (4.2) 
The purpose of the adaptive algorithm is to adjust the filter coefficients of the adaptive 
linear combiner to minimize the mean-square error. A general expression for mean 
square error as a function of the filter coefficient values, assuming that the input 
signais and the desired response are statistically stationary and that the filter 


coefficients are fixed, can be derived in the following manner [Ref 6]. The squared 


error 1S 


e*(k) = d@(k) - 2d(k)X!(k)B(k) + BE(k)x(k)x !(k)B(K) (4.3) 


Taking the expected value of both sides yields the mean square error, 


Efe*(k)] = E[d?(k)] - 2 E[d(kyx! (ky) Bek) + BE(k) E[X(K)X (ky B (4.4) 


Defining the vector P as the cross-correlation between the desired response and the 


input vector, we have 


P = E[ d(k) X(k) | (4.5) 


Similarly the input autocorrelation matrix B is defined as 


Si, 


RB = E[ X(k) X!(k)] (4.6) 


Thus the mean-square error can be expressed as 


e(k) = Ele2(k)] = E[d“(k)]}- 2B! Bck).+ B!(k) B B(k) (4.7) 


The gradient Ve(k) of the mean square error function is obtained by differentiating 
Eq.(4.7) with respect to the filter coefficient vector as follows: 


6E[e*(k)] 
Bb (k) 
Ve(k) = : = -2P +28 B(k) (4.8) 
GE[e7(k)] 
Ae 
ONIN | 


The LMS algorithm is an implementation of the method of steepest descent. 
According to this method, the next filter coefficient vector is equal to the present filter 


coefficient vector B(k) plus a change proportional to negative of the gradient, Ve(k): 


B(k+1) = B(k) - wp Ve(k) (4.9) 


The parameter pt is the factor that controls stability and rate of convergence. In other 
words, the first term on the right hand side consists of the past information and the 
second term represents the new, or updated information. 

The LMS algorithm estimates an instantaneous gradient in a crude but efficient 
manner bv assuming that 2*/ k), the square of a single error sample, is am estimaterom 
the mean-square error and by differentiating e*(k) with respect to B(k). ihe estimaved 


gradient is given by the following expression 
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Ge( k) Ge(k) 








Gbo(k) dbo(k) 
7e(k) = | = 2 e(k) (4.10) 
ae2(k) Ge(k) 
| abn. 1(4)| | ab. (K) 


The estimated gradient components are related to the partial derivatives of the 
instantaneous error with respect to the filter coefficient components. Thus the 


expression for the gradient estimate can be simplified to 


Ve(k) = -2 e(k) X(k) (4.11) 


Using this estimate in place of the true gradient in Eq.(4.11) yields the Widrow-Hoff 
LMS algorithm 


B(k+1) = B(k) + 2 p e(k) X(k) (4.12) 


Since the filter coefficient changes at each iteration are based on imperfect gradient 
estimates, we would expect the adaptive process to be noisy, that is, it would not 
follow the true line of steepest descent on the performance surface. The LMS algorithm 
can be implemented in a practical system without squaring, averaging, or 
differentiation and is elegant in its simplicity and efficiency. Each component of the 
gradient vector is obtained from a single data sample without perturbing the filter 


coefficient vector. 


C. ADAPTIVE LATTICE ALGORITHMS 

The parameters to update in a multistage lattice are its reflection coefficients. 
several algorithms have been proposed in the literature for updating the .attice 
reflection coefficients [Ref. 7,8,17,18,20,23-28]. In this section we briefly summarize 
some of those algorithms. Consider a lattice filter of order M. For stage m of the 


filter, the flow of signals is described by the following pair of equations. 
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Figure 4.3 A Single Stage of Lattice Filter. 


fink) isi fm- 1%) Ba Koi f 8m-1(k-!) (4.13) 


where m= 1,2, ... ,|M 


For stage m, the forward reflection coefficient 1s denoted by kone and the backward 
reflection coefficient is denoted by Kin.b' f(k) and g,,(k) are the forward prediction 
error and backward prediction error of stage m respectively. Before presenting the 
adaptive algorithms for the l!attice, let us quickly summarize some non-adaptive 
reflection coefficient estimation methods. Later on we can obtain the adaptive 


updating equations as approximations to these methods. 
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1. Non-adaptive Methods 
When the lattice coefficients have fixed, non-adaptive values, several methods 
have been proposed for computing these values as functions of the correlation statistics 
of the input. Two of these which are based on mean-square error (MSE) minimization 


(Ref.28] are given below. 


Method 1 : In this method, we choose kf tO minimize the mean forward 
prediction error power, E(f atk) and k,, , to minimize the mean backward prediction 
error power, Elg?_,(K)]. Here we give the final equations for k,, ¢ and k,,, ,, without 


going into the details. The forward and backward reflection coefficients are given by 


Effin 108m. (k- 


Kin f = 2 

Elen (K-) me 
Co Ellin 681 I 
ao BLP 1S 


where both , f and Ln are obtained from the cross-correlations between the 
3 3 


jth 


(m- order forward and backward prediction errors normalized by respective 


prediction error powers. 


Method 2: For a single channel lattice using real data and coefficients we can, 


m,f~ <m,b 
we can now minimize either Elinor or E[g? m(K)] in order to obtain an optimum K,,. 


however, show that k a Based on the condition that Kn fo Km b> km 
However, it seems more logical to minimize the sum Elf (kK) + Ele? (4) as suggested 


in [Ref. 28}. The resulting reflection coefficient equation in its final form is given by 


7 Pre (Ke {k-1)] 
Son f - Eire a eet jeans 

Bu eee Elg ee | 
where we have normalized the cross-correlation term with respect to the sum of mean 
prediction error powers. In both of these methods, the coefficients at stage m can be 


computed independently of those following that stage. Thus, optimum values can be 


4\ 


computed successively along the structure without affecting the other coefficients. This 
phenomenon lends to the modular nature of the lattice structure. 
2. Adaptive Metheds 
An instantaneous, gradient-descent adaptive algorithm minimizes a mean- 
square error criterion can be derived for a generalized problem. An adaptive lattice 
structure is suggested by incorporating time varying coefficients Km, (8) and Km, b8) in 
the lattice and generating algorithms for the two methods described above. The 


resulting procedures are : 


Method | : Corresponding to method 1 of non-adaptive procedures metioned 


above, we can obtain the following update equations 


yl 
Mek 1) ale one —s [Ean Sasa een 
m-l, (4.16) 
k (k=l) =k i(k) + = ets) Sane 
m,od m,o coe oh) m- | m 


where {ff is an adaptive step size parameter, A is a positive weighting constant, and the 


forward power estimate at the (m-1)% stage, Ce dk), iS 


6° an 1 Ch) =} 67.1 KD) + (1-2) [97 10k] 


and the backward power estimate at the (m-1)*h stage, Ca 5(K), 18 


O71 gk) = ® Fmetye(Kl) + (1-0) TF 70) 


Method 2: Eq.(4.15) can be approximated to obtain a recursive update equation 


aS Lollaws: 


Kin(k +1) = Ky(k) + Ht ff) gp) + fy 8m (K)] (4.17) 
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For the basic structure of the single-channel lattice, ie, Koreas Kon tne reflection 
coefficients or filter coefficients k,(k) may be updated using a modification of the 
LMS algorithm of the LMS algorithm [Ref. 17,20]. 


k (kt 1) = k(k) + = 


m-1‘*) Um) €m- 10-0 + fn. 100800] (4.18) 


pth 


2 . . - : 
where 07 ,_;(k) is the power estimate at the (m- stage. Now the updated power 


estimate is 


a7 (kh) = Aa? (kel) + (I-A) (p27, (kel) + F(R) (4.19) 


where dX is a positive weighting constant satisfying the criterion 0< <1 then controls 
the bandwidth of the filter and the resulting power averaging time. A power estimate 
is required at each stage in the lattice due to the fact that the forward and reverse error 
sequences have decreased power with increased stage number. 

Method 3: A third successful method of implementing an adaptive algorithm for 
FIR lattice structures has been reported by Griffiths [Ref. 17,18]. This algorithm has 
been originally discussed for a noise canceller application. The form of the lattice 
noise-cancelling adaptive filter is shown in Figure 4.4. The lattice noise-cancelling 
adaptive filter consists of an M stage linear prediction lattice for the reference signal 
XAK) together with a set of tap coefficients V_.(k) which provide the noise-cancelling 
subtraction paths. Griffiths’ algorithm is briefly presented in the following. The 


update equations for k_,(k) are given by 


ul { 
ae eS 


where the power estimate at the (m-1)%h 


stage, Ga 1 (k) 1s 


ork) = AO (Kel) + (I-A) [2 (KL) + Py C0] (4.21) 
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Figure 4.4 Lattice Form Implementation of Noise-cancelling Filter. 


Each coefficient V_(k) can be determined independently of V.(k) for n> m, because of 


orthogonalization provided by the lattice. Thus the resulting algorithm is 


ul 


Vink +1) = Vin(k) + oo 


(a(k) San(K)] (4.22) 


where associated power measurement Y? lk) 1s 


Yk) = RY (KL) + (1A) [Ba( KL) Sen) (4.23) 


h 


and €_.(k) is the m"" stage error signal as shown in Figure 4.4. 


te 


D. SIMULATION RESULTS 
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Figure 4.5 System Identification Experiment. 


For the purpose of computer simulation, let us consider an approximate 
algorithm given by 


K(k +1) = kk) + mT [Sy (K-1)] e(k) (4.24) 
m1 
where 
o? x(k) = Xo, (k-1) + (1-2) g* (KL) (4.25) 


The configuration used in the simulation is a system identification experiment as shown 
in Figure 4.5. The fixed filter transfer function considered is given by 
H(z) = 1 - 0.89 z! +025 2 

The convergence performance of the LMS algorithin (as given by Eq.(4.24)) can be 
observed by plotting the error, e(k). versus the update iteration, k. called learning 
curves. The input x(k) to doth fixed and adaptive filters is a white noise sequence with 
zero mean and unit variance. Figures 4.6 to 4.8 show the learning curves for the above 
example where we have used three different values for the adaptation constant, ph. In 
the next chapter, we derive a new algorithm for updating the reflection coefficients, 
based on the least mean square principle and the steepest descent algorithm. 


Improvement in convergence speed will be shown using the new algorithm. 
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Figure 4.6 Learning Curve (jt = 0.01). 
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Figure 4.7 Learning Curve (jt = 0.1). 
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Figure 4.8 Learning Curve (u = 0.5). 
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VY. ADAPTIVE LINEAR PHASE LATTICE ALGORITHM 


A. INTRODUCTION 

In Chapter III we dealt with the realization of fixed coefficient FIR lattice filter with 
both linear and nonlinear phase characteristics. We reviewed the basics of LMS 
algorithm and some adaptive lattice realizations reported in the literature in the 
previous chapter. The three adaptive lattice algorithms discussed are direct 
approximations of their non-adaptive counterparts. None of them estimates a gradient 
as required by the LMS algorithm. I[n this chapter we present a derivation for a new 
LMS based FIR adaptive lattice algorithm and extend it to the linear phase case as 
well. The new algorithm is then used in the estimation of spectral lines in white noise. 


Results of simulation are included. 


B. LMS ALGORITHM FOR THE FIR LATTICE 
From Eqs.(3.40) and (3.41) the FIR lattice equations can be written as 


CUR) Faestan tts eke ceo Bait ks) (S28) 
g tk) = Sm Benen ate Se ees) (5.2) 
where m= 1,2,...,M, Seal for m*#M and kin are the lattice reflection coefficients. A 


realization of Eqs.(5.1) and (5.2) is shown in Figure 5.1. The lattice input is 
X(K) = go(k)=fo(k). The output of the filter is 


y(k) = Fyg(k) a 


Therefore, ihe error emis pivensoy 


e(k) 


d(k) - y(k) (5.4) 
= d(k) - syq fyg_1(k) - Kyg Sng_y(k-D) a 


48 
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Figure 5.1 Adaptive FIR Lattice Filter. 


where d(k) is the desircd signal. The objective is to minimize the mean square error 


J = E { e* (k)} (5.5) 


The gradient of J with respect to k,, and s_, respectively, are given by 


G» 
Ca 
CoD 
ae 
~~ 
Pe 


, (3.0) 








Ci 
yo 
Cu 
Pr 


and 
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oJ 
OSy4 





= 2 Eck ie, (323) 


Then the LMS algorithm can be formulated as follows: 


| oJ 
kK(k+1) = Ce 


= K(k) aa Efe(k) ¢ 20) )} (373) 
K ok 
Pe | dy(k) 
= K(k) + 2 By e(k) ( “a, 
and 
oJ 
Sy(kt 1) = Syg(k) - Wot Fas ) 
SM 
syg(k) + 2B, Efe(k) fy _1(W)} (5.9) 


rr 


Where Hy, and H, are the adaptation constant, and we have replaced the true gradient 


by its instantaneous estimate. Defining 





2k) = (SE 


vields 


k(k+1) = K(k) + 2 By e(k) 2(k) 


(Sum 
sy(kKt1) = syg(k) + 2 He e(k) fy (K) 
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Where j=1,2,....M. The next step is to estimate the gradient vector z(k). For this 
purpose, define 





y = 2 512 
ij) = ok, (9.12) 
and 
dg-(k) 
F; (Kk) = 3k (bs) 
Then from Eq.(5.3) z(k) is given by 
2) = Pyj i 


(5.14) 
where 944 7 (Ck), /Ok;). Substituting Eqs.(5.1) and (5.2) to Eqs.(5.12) and (5.13) then 
we have 

®; (k) = s(k)- me ~ 


Og:_)(k- 1) 

4 ar +g eng 
ok; 8}- 

er ke I k 
© = so — a + kik) Bi + 10g 


therefore. 


®; i(k) = s:(k)@;_ 1 jf) + ki(k) Ps j ‘(k-1l) + gi. (K-1)9; | (5.15) 
Fi j(k) = si(k) 0; -1,j& ee kK (k)®;_ 1) + f- -1(K)9; ; 
where 


(5.16) 


S| 


6; | om Bk. (S307) 


is the Kronecker delta function, and i=1,2, ... ,M and j=1,2, ... ,M. If 1=0, then 
Eqs.{5.15) and (5.16) are 


© — Ox{k) | 
> re (5.18) 
J 


where j= 1,2, ... ,M. Figure 5.2 shows the computation of gradient elements for the 
adaptive FIR lattice algorithm. 


From the Eqs.(5.15), (5.16) and (5.19) we have 


D700) = Yee (5.20) 


where 1SiSj-1. The computations of gradient elements at each case of ;}= M,M-l, ... 
1 are as follows: 


Case j=M: 
OMe Sonics eames 
eS) = SHRIOL  yCR) + KCK) yg(K-L) + gy (R- 10; v9 


Pik) = si(k) Pid v(k-D te Ki(K)®7_ 1 i (’) oe fi. 100i v4 


where i= 1,2, ...,M. From Eq.(5.20), every gradient element equals to zero except final 
stage elements and gradient elements of the final stage are 
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Figure 5.2 Computation of Gradient Elements for the Adaptive FIR Lattice Algorithm. 


Ove MOS) = Sy(W@ yey w(K) + KyiK yy MiK LD) + @yy-y(k-DOyy vy 
PM(k) = Sy (kPa MRD + Ky (KO yp mC) + fy. (Ou os 


Applying the Eqs.(5.17) and (5.20) yields 


a. 


Dy Mlk) = 8y-y(k-1) (5.21) 
Fumo) = fyi) (5.22) 


Case j= vier 
Pou.) = Foys(k) = 0 
®i 1K) = S)Oi 7 1G) + KP yep) + 35K DO; yyy 
Pi Mik) = Si) Pip yD + KG )P iy yg) + 6.5098 yy 


where 1= 1,2, ... ,M. Applying the Eq.(5.20), last two stage elements are considerable, 
if 1= M-1, then we have 


Dyy-1M-1(8) = S8yy.p()Pyy.2 4.1K) + kay. (k) PF yg-2 wy D + 8yy.olk-DOy 1 Mey 


Pv-1.M-10) = 8y.1)'P yo yeh D + Ky @Oyy.2 v1 + fy.2)0y1 My 


Applying the Eqs.(5.17) and (5.20) yields 


®u-iM-1) = sy-2k-D Cm 


And the last stage terms are 
viM-1G9) = sylOP yy. a-10) * kKyGOF tor) 7 8y-1-Dey vt 
PyM-108) = sy) Pep Mr D) + Ky GP yp (9) + fy-(OM M-1 
Using the Eqs.(5.17), (5.23) and (5.24) yields 
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Finally, 
Case j=l: 


Dp i(k) = Yo y(k) = 0 


; y(k) = si()O}_ 1 1(K) + k(R)}y (KL) + 93.4 (K-1)8; (5.27) 
Ws Cs) = si(k)'Piy pOL) + KK)O} | (Kk) + £8; | (5.28) 


where 1= 1,2, ... ,M. 
The LMS algorithm derived in the foregoing can be summarized in Table 1. 


Simulation Results : 

The performance of the algorithm summanzed in Table 1 has been observed by 
computer simulation. The configuration used is the system identification experiment as 
discussed in Section (IV. D) using the same fixed coefficient filter. The learning curves 
obtained using the new algorithm are shown in Figures 5.3 to 5.5. Comparing these 
with the learning curves in Figures 4.6 to 4.8, which are obtained using approximate 


algorithms, we observe significantly faster convergence rate for the new algorithm. 


a5 


TABLE 1 
LMS ALGORITHM FOR THE FIR LATTICE 





aC) = <a emo ae 
= s(K)fmai(K) + K(k) dp-q( K-21) 
= s(k)dn-y(K-1) + Kp () £40k) 
m= 1, °2, «cee 

y(k) = fy(%) 


Update Equations: 
eee Cy / O° CK) e(k) z5(k) 


Sy( K+1) = sy(k) + 2 {./o7,(k)] e(k) fy_4(k) 


Hh 
oe 
a 


Q 
= 
= 
| 


where wt, and HP, are the adaptation constants, 
Og () = 6°, (Ke1) + (1-2) Oy 5(X) 
and 
o7_(k) = kh o7(k-1) + (1-2) £2y_4(K) 
Ss S M=-1 


are estimations of power in Z.(k) and fy_j(k), 
respectively and i is a positive weighting constant, 
O<) ie 


Gradient Vector Elements: 


@, 5(K) = 8y(K)@y_y (kK) + KCK) PG, (R21) 


oc one 
~ + 7s 


ad 


Tig) = $4( 5) P5-4 568-4) i 45(5)O324 50%) 
+ £521(%)0; 5 
z3(k) = ®y 3() 
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Figure 5.4 Learning Curve (wt = 0.3). 
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Figure 5.5 Learning Gurve (py — 


C. LINEAR PHASE FIR LATTICE ALGORITHM 
We now extend the LMS algorithm derived in the previous section to the linear 

phase FIR lattice filter. 
From Eq.(3.29), output of the linear phase FIR lattice can be written as 

y(k) = fyg(k) + gyg(k-D) (5.29) 
The lattice input is x(k) = fo(k) = go(k). Substituting Eqs.(5.1) and (5.2) in the output 
equation of the filter yields 

HED = SEMIS) © Kye BMD © SM 8M 1(k-D-1) 7 Kt IM-1-P}, 40) 


— $M U1) + 8y-i(k-D-1y) 7 KM U.16-D) + 8yp.1 (kD 


and the error, e(k), is given by 
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Figure 5.6 Adaptive Linear Phase FIR Lattice F liter. 


e(k) = d(k) - y(k) 


5.31) 
= dC) - syq [yy -1(®) + 84-1¢k-D-1)] > kM Ein (&-D) + evi 


where d(k) is the desired signal. A schematic of the linear phase FIR lattice realization 


is shown in Figure 5.6. From Eq.(5.31) we see that y(k) is a function of both sy, and 


Ky. But fyy_y(k) and gaq_)(k-1) are functions of kyy_) and so on. Therefore, in order 
to minimize the mean square error, we define a cost function 


ere { e- (k)} 


and a gradient element 


oo 


Oy(k) 
1 (5.33) 
Ok; 


Then following the treatment in Eqs.(5.6)-(5.11) the LMS algorithm for the linear 


phase lattice can be written as 


k(k+1) = K(k) + 2 My e(k) 2(k) 


| 5,34) 
s\g(kKt1) = syg(k) + 2 Wy e(k) Efyy.1(K) + gy] (K-D-DI . 


where }=1,2,...M. From Eqs.(5.29) and (5.30), the gradient vector z(k) is given by 


Z{k) = Ov (0) ot By j(k-D) 
= Syg(k) [Oy jf) + Peg j(k-D-D] + Kyg(k) [P yyy K-D) 
+ By eD) + (Fy gD) ~ gyp.10-DI By; 


(5.35) 


where 64 , = (eky4/¢k;), and On and F m,j are obtained on the same lines as in 
Eqs.(5.15)-(5.17). The resulting LMS algorithm is summarized in Table 2. 


Simulation Results : 

The performance of the algorithm summarized in Table 2 has been observed by 
computer simulation. The configuration used is the system identification experiment as 
discussed in Section (IV. D). We have used two different fixed coefficient linear phase 
FIR filter examples, one with N even and the other with N odd. They are 

H(z) = 0.154 + 0.462z"! + 0.46227 + 0.15427 
and 
H4(z) = 0.15 - 0.452! + 0.36272 - 0.4523 + 0.1524 


The learning curves obtained using the new algorithm are shown in Figures 5.7 to 


Css 


10 
Figures 3.7 and 3.8 are obtained iearning curves with linear phase ik ramsie. 
function H(z). Figures 5.9 and 5.10 are obtained learning curves with linear phase 
FIR transfer function H(z). 


60 


TABLE 2 


LINEAR PHASE FIR LATTICE ALGORITHM 





meee nalization: 
K(O) = 0 
Sy( 9) = 1 
Lattice: 
x(k) = fo(k) = go(k) 


fal K) Smt K)fy-30%) + KCK) 9,40 K-11) 

Onl) = SC K) 9,70 K-1) + KCK) £y_40%) 

Deeley... , M 

y(K) = sy 0 £y.40%) + Syiy(K-D-1)] + Ky [ fy_i(k-D) 
+ Gyiz(k-1)] 


Update Equations: 
kj(k+1) = ky(k) + 2 [my /o%, ()] e(k) 25() 
2 


2 
sy( k+l) = sy(k) + 2 (u,/or,(k)] e(k) [ £y-1() 


+ Gy ( k-D-1)] 
where wy, and ws are the adaptation constants, 
Z = ie Z 
o Car =X 6 eae) + (1-)) [®y 5(*) +E yy, 4( RD) I 
and 
2 a 2 7 _ -D-1)14 
Sm 1) XN 6 ea) (al d) eee) + Gy. 1(k-D 1)] 


are estimations of power in z.(k) and 

[ f£y-4()+9y.1(k-D-1)], respectively and XW is 

a positive weighting constant, O<AS1. 
Gradient Vector Elements: 


@, .(k) = $y(%) Oy 3() + RO) Hy RD) 
ie oe are 


+ £4-1()05 


J 
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Figure 5.8 Learning Curve (p = 0.3). 
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Figure 5.10 Learning Curve (u = 0.1). 
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D. SPECTRAL ESTIMATION 

In this section, we extend the linear phase FIR adaptive lattice algorithm derived 
in the previous section to the estimation of spectral lines in white noise. The spectral 
estimation problem is some what different from the system identification experiment 
that we have been considering so far. In a typical spectral estimation problem. we are 
given a time series and we need to estimate its spectrum. A suitable configuration that 
is frequently used for this purpose is shown in Figure 5.11. Comparing this with 
Figure 4.5, we notice that we have the given time series at the input of the adaptive 
filter and the filter updates its coefficients to minimize the mean square output. In 
doing this, we assume that the input process x(k) has been generated by passing a 
white noise sequence through a filter, say I{z), and the adaptive filter attempts to 
realize an inverse filter, say H(z)=1/I(z). Considering that the adaptive filter has 
sufficient degrees of freedom, the output of the filter will be a white noise sequence. 

The adaptive algorithm summarized in Table 2 can be used in spectral estimation 


after appropriate modification. In the present case we minimize the cost function, 
Ee 2 
Sea 


rather than J = Efe(k)}. The resulting update equations can be shown to be 


k(k+ 1) = kj(k)- 2 (By/O°4 0) e(k) z;(k) (5.36) 
sy(kt 1) = syg(k) - 2 [,/6,(k)] e(k) [yg p(k) + gy. y(K-D-1)] (5.37) 


where Z: i(k), Be (k) and o(k) are as defined in Table 2. 
Using the adaptive algorithm in Eqs.(5.36) and (5.37) we can obtain values of the 
reflection coefficients and the s,, coefficient. We then obtain the equivalent 


polynomial, B(z), bv going through the standard step up procedure given by [Ref. 1]: 


0 


ma, | © 
bam = Km 
Omi 2 Omi F Ky Om-1,m-i 
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Where 1=1,2,...m-l and m=1,2,....M. Finally the required linear phase transfer 


function is obtained as follows: 
Hy (2) = Fa(z) + 2°? Gyi(2) 


where F(z) = B(z) and G2) zMByz}), The spectrum is computed as follows: 


s) = —— 
[hp thy eH + + hay, ecHar DO 2 


where @=2nf, and f is the normalized frequency (with respect to the sampling 
frequency) in the range, OS f<0.5. 


Simulation Results : 
The input process x(k} consists of a signal in noise, given by 
x(k) = s(k) + w(k) 
where s(k) may be a single, or multiple sinusoids and w(k) is a zero mean unit variance 
white noise sequence. In the following we consider several examples using parameters 
ranging from a single sinusoid to 4 sinusoids, SNRs from 30dB to 10dB, and filter 
order, M, from 2 to 30. 


Example 1: We consider a single sinusoid given by 


x(k) = ./2 cos(2nfk) + w(k) (5.38) 


where f=0.15, SNR = 30dB. 

mecress.12, 2.05, 5.14, and 5.15 are plots of frequency spectrum with different 
order of jattice Vi and adaptation constant Hp. 

Figure 5.12 is the plot of the case M=2 and pt =0.015. We see that the peak is 
at f=0.15 and no spurious responses. 

Figure 5.13 shows the frequency spectrum of M=4 and jt=0.01. We observe 


that one peak is at f=0.15 and a spurious response whose magnitude is about -53dB at 
f= 0.37. 
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Figure 5.11 Spectral Estumation Mode. 


Figure 5.14 shows the frequency spectrum of M=10 and }t=0.03. There are one 
peak at f=0.15 and four spurious responses at f=0.05, 0.25, 0.35 and 0.45. The largest 
spurious response is about -15 dB at f=0.45. 

Figure 5.15 shows the frequency spectrum of M=20 and jt=0.03. There are one 
peak at f=0.15 and nine spurious responses. The largest spurious response is about 
-27 dB at f=0.125, 

Through the Figures 5.12, 5.13, 5.14, and 5.15 we can determine that the best 
model order to detect the one sinusoidal signal is VI = 2. 

Next, we fix the order of lattice M, adaptation constant tt and iteration number 
k. Figures 5.16, 5.17 and 5.18 are the plots of frequency spectrum with changing the 
signal to noise ratio (SNR). 

Figure 5.16 is the plot of M=20, n=0.05, k= 3000 and SNR=30dB. There is a 
peak at f=0.15 and the largest spurious response is about -27dB at f=0.125. 


66 


Figure 5.17 is the plot of M= 20, w#=0.05, k= 3000 and SNR=20dB. The largest 
spurious response at f=0.425 is larger than the desired response. At this case, we 
cannot estimate the frequency of sinusoid. 

Figure 5.18 is the plot of M=20, #=0.05, k=3000 and SNR=10dB. Three of 
the spurious responses are larger than the desired response. The desired response is 
about -13dB. 

From Figures 5.16, 5.17 and 5.18, if the SNR is getting smaller then estimating 


sinusoidal frequency is more difficult. 


Example 2: Next, to estimate two sinusoidal frequencies in noise, set the input x(k) as 


x(k) = /2 { cos(2mfyk) + cos(2mf>k) } + w(k) (5.39) 


where the normalized frequencies of signals fj =0.15 and f,=0.25 and set SNR = 30dB. 

Figure 5.19 shows the frequency spectrum of M=4 and t=0.02. There are two 
peaks at f; =0.15 and f,=0.25 and no spurious responses. 

Figure 5.20 shows the frequency spectrum of M=20 and =0.022. There are 
two peaks at f) =0.15 and f,=0.25, and eight spurious responses. The largest spurious 
response is about -20dB at f= 0.375. 

From Figures 5.19 and 5.20, the best model order to estimate the frequencies of 


two sinusoidal signals is M=4. 


Example 3: Next, to estimate four sinusoidal frequencies in noise, set the input x(k) as 


x(k) = 4/2 { cos(2mf,k) + cos(2mf>k) + cos(2mf,k) + cos(2mfyk) } + w(k) (5.40) 


where the normalized frequencies of signals {, = 0.05, [,=0.15, f,=0.25, [,=9.35, and 
set SNR = 30dB. 

Figure 5.21 shows the frequency spectrum of M=8 and #=0.015. There are four 
peaks at f; =0.07, f5 =0.15, f3 =0.27 and fg=0.375, and no spurious responses. 

Figure 5.22 shows the frequency spectrum of M=30 and n=0.014. There are 
four peaks at f; =0.05, f,=0.15, f3=0.25 and fy=0.35, and eleven spurious responses. 
The largest spurious response is about -23dB at f=0.45. 
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From Figures 5.21 and 5.22, the best model order to estimate the frequencies of 


four sinusoidal signals is M = 8. 
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Figure 5.12 Frequency Spectrum (I sinusoid, M=2, p#=0.015). 
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Figure 5.14 Frequency Spectrum (1 sinusoid, M= 10, p = 0.03). 
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Figure 5.15 Frequency Spectrum (1 sinusoid, M = 20, 2 =0.03). 
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Figure 5.16 Frequency Spectrum (1 sinusoid, M=20, SNR= 30dB). 
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Figure 5.18 Frequency Spectrum (1 sinusoid, M=20, SNR= 10dB). 
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Figure 5.19 Frequency Spectrum (2 sinusoids, M =4, 2, =0.02). 
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Figure 5.20 Frequency Spectrum (2 sinusoids, M = 20, p= 0.022). 
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Figure 5.21 Frequency Spectrum (4 sinusoids, M=8, #=0.015). 
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Figure 5.22 Frequency Spectrum (4 sinusoids, M = 30, p=0.014). 
73 


E. SUMMARY, CONCLUSION AND SUGGESTED FUTURE WORK 

The thesis investigates the application of lattice structures in Prony’s method of 
spectral line estimation. The complete solution to Prony’s method consists of three 
steps: (1) representing a given process of M sinusoids in noise in terms of complex 
exponentials, (11) finding roots of a symmetric polynomial, and {i1) estimating the 
frequency, phase and amplitude information. In this thesis, however, we have 
emphasized only the frequency estimation problem. 

The underlying theory of Prony’s method has been briefly reviewed in Chapter II. 
By using a modified Prony’s approach, we observed that the frequency estimation (as 
part of Prony’s method) requires a polynomial with a linear phase property. The 
basics of the linear phase FIR filter and the lattice structure have been included in 
Chapter II]. Also we have shown that a linear phase FIR transfer function can be 
realized from an FIR lattice using D number of additional unit delays and an adder, 
where D is determined by the order of the linear phase polynomial. 

The principles of adaptive filtering and the LMS algorithm have been briefly 
studied in Chapter [V. We have derived a new LMS algorithm for the FIR non-linear 
phase ans linear phase transfer functions in Chapter V. The problem spectral 
estimation using the new adaptive algorithm has been addressed, and the simulation 
results are included. 

The significant outcome of the work is the derivation of an LMS type algorithm 
for a linear phase lattice structure and its application to spectral line estimation as part 
of Prony’s method. As has been shown, the new algorithm yielded faster convergence 
rate compared to previously reported approximate algorithms. The application of this 
algorithm to spectral estimation produced high spectral resolution as illustrated by the 
simulation results. 

However, we have observed spurious spectral responses, when the model order is 
higher that the required. Also we have observed that the SNR performance of the 
algorithm needs to be improved. For input SNRs less than about 10 dB, the algorithm 
performed poorly. Finally, we have dealt with only the frequency estimation part of 
Prony’s method. [n a future work, some of the above mentioned problems may be 


taken up. 
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APPENDIX A 


EXAMPLES OF OBTAINING A LATTICE STRUCTURE FOR THE FIR 
TRANSFER FONCTION 


Example A.1 
Obtaining a lattice structure for the general FIR transfer function. The unit 
sample response of a FIR transfer function is given by 
H(z) = 0.5 + 0.25274 + 0.1252% + 0.062529 


Solution: Given unit sample response is 
H3(z) = 0.5 + 0.2527! + 0.125272 + 0.062523 Can 
Using Eq.(3.22), we have polynomial for 3 sections 
H3(2) = 03,9 + 5,27) +b3,.92°° + 03 32° od) 


Comparing Eqs.(A.1) and (A.2), we have 


b3 9 = 0.5 
b3 4 = 0.25 
(A.3) 
Starting with m=3 we have from Eq.(3.46) 
(A.4) 


$3 = Sages 0.5 


Now, we need to generate the coefficients for H(z) and from Eq.(3.46) 
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b i SP m.i Km m-i 
m-1],1 S 2 - k 2 
m m 


And for m=3 and i=0 we have 


q 
$303 (1 + K3b3 3 


oe 
_ _ (0.5)(0.5) - (0.0625)(0.0625) _ 
2,0 (0.5)* = (0.0625) 


Next, form = 3 andi = | we get 


ee esSe salmaaa 
om 
| $3° * K3" 


Oro (0:25 ea OGes (O25 
by, = en Sahel =" @gens 
(0.5)* - (0.0625) 
and finally 


_ $363 9 - k3b3 | 


b 
a S35 ee 
0.5)(0.125) - (0.0625)(0.25 
sete (0.5) , eat X ) = 9.19048 
(0.5)* - (0.0625) 


from Eq.(3.46), we have 


0.19048 


YN 
o) 
I 
or 
td 
to 
iI 


The new polynomial is 


ee -| -2 
H4(Z) = b 0 te bo 42 —- ba 22 
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(A.5) 


(A.6) 


(A.7) 


(A.8) 


(A.9) 


H4(z) = 1 + 0.476192"! + 0.190482" 


for m= 2 and i=0, we have 


mee 220" 822.2 
L. 


{ - Ks 


; 1 - (0.19048)(0.19048) 
L,I 1 - (0.19048) 


4 
& 


and finally 


Bo 80> | 
eo 


0.47619 - (0.19048)(0.47619) 
b) |) = ——— . = (2.40000 
! L - (0.19048)- 


From Eq.(3.46), we have 


Ss] = d10 = 1 


(A.10) 


(A.11) 


(A.12) 


(A.13) 


Therefore, reflection coefficients and s coefficients of the FIR lattice structure are 


k) = 0.40000 
ky = 0.19048 
k; = 0.0625 
Sy = | 

39 L 

33 = 0.5 
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Example A.2 (N= 4:even) 


Obtaining a lattice structure for the linear phase FIR transfer function. A linear 
phase FIR transfer function is 


H4(z) = 0.154 + 0.46227! + 0.462272 + 0.15429 (A.14) 
Solution: Eq.(A.14) can be written in the form of Eq.(3.35) 
7 = -1 a 2-2 -! AiG 
H3(Z) = a9 faye Zim cua yz) (A.15) 
Comparing Eqs.{(A.14) and (A.15), we get the following relationships. 
ao) = 0.154 
(A.16) 
ay = 0.462 


To determine the lattice reflection coefficients of the corresponding linear phase FIR 
lattice structure and the number of unit delay, we use Eq.(3.29) which 1s 


Hy. (2) = Fyy(z) + 2? Gyy(z) (A.17) 


From Eq.(A.14), we have N=4(even), the order of the polynomials Fy4(z) and Gyy(Z), 
M = (N/2)-1 = (4/2)-1 = 1, and the number of unit delays, D = N/2 = 2. Now, 


from Eq.(3.37), we may have the forward prediction error transfer function as follows: 
F(Z) a ay te ay zh 
(A.18) 
= 0.154 + 0.462 z’! 
[Om Vi= 1, 2d4 A218) can be remmiten es 
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From Eqs.(A.18) and (A.19), we get the desired reflection coefficient and s coefficient 
of the linear phase FIR lattice structure as follows: 
Ky = bi = (0.462 


Sj a big = 0.154 
Example A.3 (N= 5:0dd) 


Obtaining a lattice structure for the linear phase FIR transfer function. A linear 
phase FIR transfer function is 


H(z) = 0.15 - 0.4527! + 0.3622 - 0.45273 + 0.1524 


(A.20) 
Solution: Eq.(A.20) can be written in the form of Eq.(3.30) 
ES -| 7 -2 -2 -1 -2 , 
H(Z) ay oe ayZ te (1/2)a5z Toa g, {(1/2)a5 a az a ayZ } (A.21) 
Comparing Eqs.(A.20) and (A.21), we get the following relationships. 
a, = -0.45 (A.22) 
a4 = 0.36 


To determine the lattice reflection coefficients of the corresponding linear phase FIR 
lattice structure and the number of unit delay, we use Eq.(3.29) which is 


Hy.) = Fry(2) + 2°? Gy(2) (A.23) 


ee) 2 


From Eq.(A.20), we have N= S(odd), the order of the polynomials Fy,4(z) and Gyj(Z), 
M = (N-1)/2 = 


2. and the number of unit delays, D = (N-1)/2 = 2. Now, 
from Eq.(3.52), we may have the forward prediction error transfer function as follows: 
F(z) = ag +a, zi +(1 ~ 
(Z) ag + ay Z (1/2) a5 Z 


(A.24) 
0.15 -0.45z! + 0.18 2 
for M=2, Eq.(A.24) can be rewritten as 


fe 


a -] -2 
F(Z) = b2 0 a b2 | Tae os b 2 Zz, 


From Eqs.(A.24) and (A.25), we have 
k, = bo 2 = 0.18 
34 = D4 0 = 0.15 


Now, we need to generate the coefficients for F(z) and from Eq.(3.46) 


b 
1,0 0.157 - 0.18 


, pee 


_(0.15)(-0.45) - (0.18)(-0.45) 


Syb4 1 - kyb 
bp, = et = -1,36364 


: aoe (GS 0.18 


From Eqs.(A.26) and (A.27), we have 
aa bY = -1.36364 
S] = dy 0 = | 


(A.25) 


(A.26) 


(A.27) 


Therefore, reflection coefficients and s coefficients of the linear phase FIR lattice 


STIUCTUTe ane 


k) = -1.36364 
Ss] = J 
9) = 0.15 
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APPENDIX B 
COMPUTER PROGRAMS 


(Ce Pes eee ee yee Ce ee) eee era haere te MB tages eee ie Oe Ppa wee he SUL les CEP ee wy MRL ow br Pe) ee. Oe 


* This program is designed for the system identification experiment, * 
* which is shown in Section (IV.D). The learning curves can be obtained* 
* by plotting the error, E, versus the update iteration, k. * 


KRAEKKKKKKKKKKKKKKAKKKKKKKKRKKRRKRRRKRRKRARKRRKKRKRRKRKRRKRRRKRRRKRRRKRRRKRRRRRKRERARKKR 


Variable Definition 


ISEED : seed for the random number generation (white noise) 


AMU : adaptation constant 

k ; time index 

M : order of the FIR transfer function or total number of 
lattice sections 

NB ; number of iterations 

A(I) : coefficients of the FIR transfer function 

B(I) : reflection coefficients of the lattice 

F(I) : forward prediction error 

G(I) : backward prediction error 


GD(I) : delayed backward prediction error 
SGL(I): estimations of power 

W(K) : input of both FIR and lattice 
YF(K) : output of the FIR filter 

YL(K) : output of the lattice filter 
ER(K) : squared error 


Variable Declaration 


+ OF OF OF OOOH HH 


PVE CERetonnD,K,1,7,N,M,NB 
REAL A(100),B(100) ,F(100) ,G(100) ,GD(100) , YF(10000) , YL(10000) 
REAL X(100),SGL(100),ER(10000) ,W(10000) ,¥(10000) ,AMU,E 
Set Adaptation Constant uw and Number of Iterations NB 
ee oo) 
i FORMAT (5X, 'AMU',5X,'NB') 
READ(5,%*) AMU,NB 


* inpera lization 


DO 10 K=1,100 
A(K)=0 


8] 


+ + + + 


10 


iS 


20 


40 


45 
30 


+ 


B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=6 
SGL(k)=1.0 
continue 
DO 15 K=1,10000 
YF(K)=0 
YL(K)=0 
ER(K)=0 
W(K)=0 
Y¥(k)=0 
CONTINUE 
E=0. 
R=1 
ISEED=343169 
M=2 
A(1) 
A(2) 
A3) 


1.0 
= One 
TORS 


Random Number Generation 
(mean zero, unit variance, white sequence) 


DO 20 N=1,NB 
CALL LNORM(ISEED,RN,1,1,0) 
W(N) = RN 
CONTINUE 
FIR Filter Calculation 


DO 30 K = 1,NB 


X(1)=W(K) 

YF(K)= A(1)*X(1) 

DO 40 I=1,M 
YF(K)=YF(K)+A(I+1)*X(I+1) 

CONTINUE 

DO 45 I=1,M 
X(M+2-I)=K (M+1-1) 

CONTINUE 

CONTINUE 


Lattice Filter Calculation 


DO 50 K =1,NB 
F(1)=W(K) 
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+ 


+ 


60 


70 


50 
300 


200 


8 


G(1)=W(K) 

DO 60 I=1,M 
F(I+1)=F(I)+B(I)*GD(T) 
G(I+1)=GD(I)+B(I)*F(I) 

CONTINUE 

YL(K)=F(M+1) 


Calculating the Error 
E = YF(K) - YL(K) 

Updating the Reflection Coefficients 
CALL LMS (B,GD,E,AMU,SGL,M) 


ER(K)=E*%*2 
DO 70) J=1 ,M 
GD(J)=G(J) 
CONTINUE 
IF (K.EQ.R) THEN 
WRITE(6,300) K, B(1),B(2),ER(K) 
R=R+50 
END IF 
Y(k)=E 
CONTINUE 
FORMAT (3X,16,5X,3(F10.7,4X) ) 


Plotting the Learning Curve 


CALL PLOT(Y,N) 
STOP 
END 


SUBROUTINE LMS(B,GD,E,AMU,SGL,M) 
REAL B(100),GD(100) ,SGL(100),E,AMU 
INTEGER 
DO 200 I=1,_ 
SGL(I)=0.9*SGL(I)+0.1*(GD(I)*GD(I)) 
CONTINUE 
DO 210 I=1,M 
B(I)=B(1I)+(AMU/SGL(I))*E*GD(I) 
CONTINUE 
RETURN 
END 


§3 


10 


SUBROUTINE PLOT(Y,N) 
DIMENSION Y(N) ,X(10000) 
ISTP=N/10 

DO 10 J=1,N 


X(J)= 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


J 

TEK618 

PRTPLT (72,6) 

SHERPA ('PARKPARK','A' ,3) 
PHYSOR(1.,1.) 

RESET('ALL') 

PAGE(8.5,11.0) 
HWROT('AUTO') 

XINTAX 

AREA2D(5.0,2.8) 

HEIGHT (0.10) 

COMPLX 
SHDCHR(90.0,1,0.002,1) 
HEADIN('LEARNING CURVES$',100,2.0,1) 
XNAME ( 'ITERATIONSS' ,100) 
YNAME ('ERROR$' ,100) 
MESSAG('ADAPTIVE LATTICE(AMU=.5,FIG4.8)$',100,3.0,-0.8) 
FRAME 
GRAF(0,ISTP,N,-3.0,1.5,3.0) 
THKCRV (0.02) 

CURVE (X,Y,N,0) 

ENDPL(0) 

DONEPL 


RETURN 


END 
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AAKAKAP rogram QRAKKKAKAAAKAKKAKKAKKAAKRAKKKRARKAAKKAKKKKRRKKARKRRKRKRRKRKRKKRKKRAK 


* This program is designed for the system identification experiment * 
* using the LMS algorithm which was derived in Section (V.B). * 
KAKRKARKAKKRKKRAKKAKARKARKKA AKA RAK KKK ARK AR ER AKA RRR AKKRAKAKA KAR RARAKAKKAKARKAR 
: : 

~ 

~ 


BNAEGER ISEED,K,1,J,N;M,NB 

REAL A(100),B(100) ,F(100) ,G(100) ,GD(100) ,x(100),YF(5000) , YL(5000) 
REAL ER(5000) ,W(5000) ,SGL(100) ,PH({100,100) ,PS(100,100) 

READY PSD(100,100) ,GR(100),Y¥(5000) ,AMU,E 


* Set Adaptation Constant [f and Number of Iterations NB 


WRITE(5,1) 
1 FORMAT(5X, 'AMU',5X,'NB') 
READ(5,*) AMU,NB 


x Miacialigation 


E=Q . 
R=1 
DO 5 K= 1,5000 
W(K)=0 
YF(K)=0 
YL(K)=0 
ER(K)=0 
Y(K)=0 
5 CONTINUE 
DO 10 K=1,100 
A(K)=0 
B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=0 
SGL(K)=1.0 
GR(K)=0 
10 CONTINUE 
DO 15 K=1,100 
DO 15 L=1,100 
PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 
15 CONTINUE 
ISEED=343169 
M=2 
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>, Nr BP, wo 


40 


45 
30 


60 


70 


A(1) 
A(2) 
A(3) 


io 
-0.89 
+0.25 


Random Number Generation 


(mean zero, unit variance, white sequence) 


DO 20 ‘N=1,NB 


CALL LNORM(ISEED,RN,1,1,0) 


W(N) = RN 
CONTINUE 
FIR Filter 


DO 30 K = 1,NB 


X(1)=W(K) 
YF(K)= A(1)*X(1) 
DO 40 I=1,M 
YF (K)=YF(K)+A(I+1)*X(I+1) 
CONTINUE 
DO 45 I=1,M 
KX (M+2-I)=K(M+1-1) 
CONTINUE 
CONTINUE 


DO 


Lattice Filter 


50 K =1,NB 

F(1)=W(K) 

G(1)=W(K) 

DO 60 I=1,M 
F(I+1)=F(I)+B(I)*GD(T) 
G(I+1)=GD(I)+B(1I)*F(I) 

CONTINUE 

YL(K)=F (M+1) 

E = YF(K) - YL(K) 


Updating the Reflection Coefficients 


CALL MLMS (PH,2S,P?SD,GR,F,G,GD,8,SGL,£,AMU,M) 
ER(K)=E**2 
DO 70 J=1,M 
GD(J)=G(J) 
CONTINUE 
IF (K.EQ.R) THEN 
WRITE(6,300) K, B(1),B(2),B(3),B(4),B(5),B(6),B(7),B(8),E 
R=R+30 
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ao 
300 
100 


+ 


+ 


+ 


10 


20 
200 


210 


21 


END IF 
Y(K)=E 
WRITE(6,100) K,W(K),YF(K),YL(K),E,ER(K) 
CONTINUE 
FORMAT (1X,16,2X,9(F10.7,2X)) 
FORMAT (3X,13,2X,5(F10.7,2X)) 


Plotting the Learning Curve 


CALL PLOT(Y,N) 
STOP 
END 


SUBROUTINE MLMS(PH,PS,PSD,GR,F,B,BD,R,SGL,EK, AMUN) 
DIMENSION PH(100,100) ,PS(100,100) ,PSD(100,100) ,F(100) ,B(100) 
1,BD(100) ,R(100),GR(100),SGL(101) 
GR(N)=BD(N) 

DO 200 I=2,N 

PH(I,1)=BD(N+1-1) 

PS(I,1)=F(N+1-I) 

ia=1-1 

DO 10 K=1,IT 
PH(I,K+1)=PH(I,K)+R(N+1-I+K)*PSD(I,K) 
PS(I,K+t1)=PSD(I,K)+R(N+1-I+K) *PH(I,K) 
DO 20 K=1,IT 

PSD(I,K)=PS(I,K) 

CONTINUE 

hor210 K=2,N 

GR(N+1-K)=PH(K,K) 

DO 211 K=1,N 
SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.0 
DO 220 I=1,N 
R(I)=R(1)+(AMU/SGL(I) ) *EK*GR(I) 
IF(R(I).GE.1.0) R(I)=0. 

PAR eee =d.0) R(I)=0. 

CONTINUE 

RETURN 

ZND 


SUBROUTINE PLOT(Y,N) 
DIMENSION Y(N),X(5000) 
ISTP=N/10 

DO 10 J=1,N 
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10 


X(J)= 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


i 

TEK618 

PRIPLT (72,6) 
SHERPA('PARKPARK','A',3) 
PHYSOR(1.,1.) 

RESET('ALL') 

PAGE(8.5,11.0) 

HWROT('AUTO' ) 

XINTAX 

AREA2D(5.0,2.8) 

HEIGHT(0.10) 

COMPLX 
SHDCHR(90.0,1,0.002,1) 
HEADIN('LEARNING CURVES$',100,2.0,1) 
XNAME('! ITERATIONS$' ,100) 
YNAME(' ERRORS! ,100) 
MESSAG('ADAPTIVE LATTICE(AMU=.5,FIG5.5)$',100,3.0,-0.8) 
FRAME 
GRAF(0,ISTP,N,-2.5,1.25,2.5) 
THKCRV (0.02) 

CURVE (X,¥,N,0) 

ENDPL(0) 

DONEPL 


RETURN 


END 
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AAAAKAAPrOgram ZAKKAKAKKAAKAAKKAAKKAAKKAAKKAKKAKKAKKAKKKKAKKAKKKAKAKAKKAKAKAAKARKAKKAAKARKKK 


* This program is designed for the system identification experiment. x 
* The LMS algorithm shown in Section (V.C) was extended to the linear * 
* phase FIR lattice filter. * 


FFI IKKKKKKKRERRKAKKKKKKKKAKRKKKKRKEKKREKKRAKKKKRAKKKEKRKKAKKKKAKKRKAKKKAKKKAKKEHK 
* 
* 
* 
INTEGER ISEED,K,1I,J,N,M,NB,R,MA,ML 
REAL A(100),B(100) ,F(100) ,G(100),GD(100) , YF(3000) , YL(3000) 
REAL X(100) ,H(100) , YD(3000) ,ER(3000) ,W(3000) , YFF(3000) , YLL(3000) 
REAL SGL(100) ,PH(100,100),PS(100,100) ,PSD(100 ,100) 
REAL GR(100),GRD(100,100),GRB(100),¥(3000) ,AMU,E,C 


+ 


Set Adaptation Constant ff and Number of Iterations NB 


WRITE(5,1) 
1 FORMAT (5X, 'AMU',5X,'NB') 
READ(5,*) AMU,NB 


Gmveve lization 


E=0. 

R=1 

M=4 

ISEED=343169 

Domo. K=1, 100 
A(K)=0 
B(K)=0 
X(K)=0 
F(K)=0 
G(K)=0 
GD(K)=0 
H(K)=0 
SGL(K)=1.0 
GRU =O 
GRB(K)=0 
10 CONTINUE 
9O 11 X=1,3000 

W(K)=0 
YF(K)=0 
YL(K)=0 
YFF(K)=0 
YLL(K)=0 
YD(K)=0 
ER(K)=0 
¥(K)=0 
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ie 


i 


+ tb 


20 


A- 


Za 


Ze 


a 


CONTINUE 
DO 15 K=1,100 
DO 15 L=1,100 


PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 
GRD(K,L)=0 
CONTINUE 
C=.154 
A(1)=.154/c 
A(2)=.462/C 
=,15 
A) =e © 
A(2)=-.45/C 
A(3)=.36/C 


Random Number Generation 
(mean zero, unit variance, white sequence) 


DO 20 N=1,NB 
CALL LNORM(ISEED,RN,1,1,0) 
W(N) = RN 

CONTINUE 

IF (MOD(M,2).E0.0) MA=M/2 

IF (MOD(M,2).NE.0) MA=(M+1)/2 

IF (MOD(M,2).EQ.0) ML=MA 

IF (MOD(M,2).NE.0) ML=MA-1 


Separation of Coefficients 


IF (MOD(M,2).EQ.0) THEN 
DO 21 I=1,MA 
H(I)=A(T) 
H(MA+2+I)=A(MA+1-I) 
CONTINUE 
H(MA+1)=A(MA+1)/2 
H(MA+2)=A(MA+1)/2 
END IF 
Lf (MOD (CM 2) E20 )) FHEN 
DO 22 t= le 
H(I)=A(T) 
H(MA+I)=A(MA+1-I) 
CONTINUE 
END IF 


LINEAR PHASE FIR FIUrTER 


90 


40 


41 


42 


43 


45 
30 


a 


60 


70 


DO 


30 K = 1,NB 
X(1)=W(K) 
YF (K)=H(1)*X(1) 
IF (MOD(M,2).EQ.0) THEN 
DO 40 I=1,MA 
VF (K)=YF(K)+H(I+1)*X(I41) 
CONTINUE 
YFF(K)=YF(K) 
DO 41 I=1,MA+1 
YFF(K)=YFF(K)+H(MA+1+I )*X(MA+I) 
CONTINUE 
END IF 
IF (MOD(M,2).NE.0) THEN 
DO 42 I=1,MA-1 
YF (K)=YF(K)+H(I+1)*X(I+1) 
CONTINUE 
YFF(K)=YF(K) 
DO 43 I=1,MA 
YFF(K)=YFF(K)+H(MA+I) *X(MA+T) 


CONTINUE 
=ND IF 
DO 45 I=1,M 
X(M+2-I )=X(M+1-1) 
CONTINUE 
CONTINUE 


DO 


LATTICE FILTER 


50 K =1,NB 
F(1)=W(K) 
G(1)=W(K) 

DO 60 I=1,ML 
F(I+1)=F(I)+B(I)*GD(I) 
G(I+1)=GD(I)+B(1)*F(I) 

CONTINUE 

YL(K)=F(ML+1) 

YD (K)=G(ML+1) 

YLL(K)=VYL(K)+¥D(K-MA) 

ae) - LLL CK) 


Updating the Reflection Coefficients 


CALL MLMS (P,Q,QD,GRB,GRD,F,G,GD,B,SGL,E,AMU,ML,MA) 


ER(K)=E**2 
pO 70 J=1,ML 
GD(J)=G(J) 

CONTINUE 


9) 


+ 


50 


300 


PEO 


120 


200 


Z10 


220 


270 


Y(K)=E 
WRITE(6,300) K,Y(K) 
CONTINUE 


Plotting the Learning Curve 


CALL PLOT(Y,N) 
FORMAT(3SX, EG Sx h10 03) 
STOP 

END 


SUBROUTINE MLMS(PH,PS,PSD,GRB,GRD,F,B,BD,R,SGL,EK,AMU,N,MA) 
DIMENSION PH(100,100),PS(100,100) , PSD(100,100) ,GRD(100,100) 
1,GRB(100) ,F(100),B(100) ,BD(100) ,R(100) ,GR(100) ,SGL(101) 


GR(N)=BD(N) 

GRB(N)=F(N) 

DO 200 I=2,N 

PH(I,1)=BD(N+1-I) 
PS(I,1)=F(N+1-1) 

IT=I-~1 

DOT OeK—l ir 
PH(I,K+1)=PH(1I,K)+R(N+1-I+K) *PSD(I,K) 
PS(I,K+1)=PSD(I,K)+R(N+1-I+K)*PH(I,K) 
DO 120 K=1,1T 

PSD(1I,K)=PS(I,K) 

CONTINUE 

DO 210 K=2,N 

GR(N+1-K)=PH(K,K) 
GRB(N+1-K)=PS(K,K) 

DO 220 K=1,N 

DO 220 L=1,MA 
GRD(K,MA+2-L)=GRD(K,MA+1-L) 

DO 230 K=1,N 

GRD (K,1)=GRB(K) 

DO 240 K=1,N 

GRB(K)=GRD(K,MA+1) 

DO 250 3= 

GR(X)=GR(K)+GRB(K) 

DO 260 K=1,N 
SGL(K)=.90*SGL(K)+.10*GR(K)*GR(K)+1.0 
DO 270 I=1,N 
R(I)=R(I)+(AMU/SGL(I) ) *EK*GR(I) 
IF(R(I).GE.1.0) R(I)=0. 
IF(R(I).LE.-1.0) R(I)=0. 
CONTINUE 
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a 


10 


RETURN 


END 


SUBROUTINE PLOT(Y,N) 
DIMENSION Y(N),X(3000) 


ISTP= 


N/10 


DO 10 J=1,N 


K(J)= 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


a 

TEK618 

PRTPLT(72,6) 
SHERPA('PARKPARK','A',3) 
PHYSOR(1.,1.) 

RESET('ALL') 

PAGE(8.5,11.) 

HWROT('AUTO') 

XINTAX 

AREA2D(5.0,2.8) 
HEIGHT(0.10) 

COMPLX 
SHDCHR(90.0,1,0.002,1) 
HEADIN('LEARNING CURVES!',100,2.0,1) 
XNAME ('ITERATIONSS' ,100) 
YNAME ('ERRORS' ,100) 
MESSAG('ADAPTIVE G-M LATTICE(F5.10,AMU=0.1)$',100,3.0,-0.8) 
FRAME 
GRAF(0,ISTP,N,-8.0,4.0,8.0) 
THKCRV(0.02) 

CURVE(X,Y,N,0) 

ENDPL(0) 

DONEPL 


RETURN 


END 
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RAAKRKAPrOgr AM GRARRAARKARKARKARKK RAK KAAKKAKERRKK RK RK AK RAK ER ARK RRA AK KKKKER 


* 
* 
x 
x 


This program is designed for the estimation of spectral lines in 

white noise. The input process x(k) consists of a signal in noise, 
and a signal may be a single, or multiple sinusoids. The algorithm 
was derived in Section (V.D). 


* 
* 
* 
* 


KRRRERKEKKRRRRRKRRKRKKRRRRRRKRRERRRRKKRRRKERRRKRRRERERRRRRRRRERRRRRRKRKRRRKRRKRKRKERK 


* 
* 
* 


INTEGER ISEED,K,I,J,N,M,NB,R,MA,ML,D 

REAL A(100),B(100),F(100) ,G(100) ,GD(100) , YL(5000) 

REAL Y(100),YD(5000) ,W(5000) , YLL(5000) , YD(5000) , INP (5000) 
REAL SGL(101),PH(100,100),?S(100,100) , PSD(100,100) 

REAL GRD(100,100) ,GRB(100) ,GR(100) 

REAL RE(100),IM(100) ,AJ(100) ,SD,SNR,AVG,AMP ,AMU,E 


Set Adaptation Constant w and Number of Iterations NB 


WRITE(5,1) 
FORMAT (5X, 'AMU',5X,'NB') 
READ(5,*) AMU,NB 


Inveralization 


ISPEC=1000 

SNR=30. 

SD=10**(-(SNR/20)) 

AMP=SQRT(2.) 

AVG=0. 

E=0. 

R=1 

PI=3.141592654 

M=8 

MP1=M+1 

ISEED=343169 

IF (MOD(M,2).EQ.0) THEN 

MA=M/2 

ML=MA 

Biot 

MA=(M+1)/2 

ML=MA-1 

END IF 

DO 10 K=1,100 

A(K)=0 
B(K)=0 
F(K)=0 
G(K)=0 
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+ OF OF OF 


+ 


GD(K)=0 
SGL(K)=1.0 
GR(K)=0 
GRB(K)=0 
RE(K)=0 
IM(K)=0 
AJ(K)=0 
¥(K)=0 
10 CONTINUE 
DO 11 K=1,5000 
W(K)=0 
YL(K)=0 
YLL(K)=0 
YD(K)=0 
ER(K)=0 
INP (K)=0 
ie CONTINUE 
DO 15 K=1,100 
DO 15 L=1,100 
PH(K,L)=0 
PS(K,L)=0 
PSD(K,L)=0 
GRD(K,L)=0 
15 CONTINUE 


Random Number Generation 
(mean zero, unit variance, white sequence) 


DO 20 N=1,NB 
CALL LNORM(ISEED,RN,1,1,0) 
W(N) = SD*RN+AVG 
20 CONTINUE 


EAT REGE FILTER 


DO 50 K =1,NB 
AK=K-1 
INP (K)=AMP*COS (2*PI*,15*AK)+W(K) 
C INP (K)=AMP* (COS (2*PI*™. 1S"AK)+COS(2*PI*.25*AK) )+W(K) 
INP (K)=AMP* (COS (2*PI*.O5*AK)+COS (2*PI*.15*AK)+COS (2*PI*.25*AK 
*)+COS (2*PI*. 35*AK) )+W(K) 
F(1)=INP(K) 
G(1)=INP(K) 

DO 60 I =1,ML 
F(I+1)=F(1I)+B(I)*GD(I) 
G(I+1)=GD(I)+B(I)*F(I) 

60 CONTINUE 
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+ Ft FF + 


70 


80 


81 


50 


300 
301 


YL(K)=F (ML+1) 

YD(K)=G(ML+1 ) 

YLL(K)=YL(K)+YD(K-MA) 
P= Yrnn) 


Uodating the Reflection Coefficients 


CALL MLMS (PH,PS,PSD,GR,GRB,GRD,F,G,GD,B,SGL,£,AMU ME Me) 


ER(K)=E**2 
DO 70 J=1,ML 
GD(J)=G(J) 

CONTINUE 


IF (K.NE.ISPEC) GO TO 50 
IF (K.EQ.NB) THEN 

WRITE(6,300) K,INP(K),E 

WRITE(6,301) K,B(1),B(2),B(3),B(4),B(5),B(6),B(7) ,B(8) ,B(9) 
R=R+100 


Determine the FIR Coefficients from the Lattice Reflection 
Coefficients 


CALL STEPUP (A,B,ML) 
IF (MOD(M,2).EQ.0) THEN 
DO 80 I=1,M/2 
A(M+2-I)=A(T) 
A(M/2+1)=2*A(M/2+1) 
ELSE 
DO 81 I=1,(M+1)/2 
A(M+2-I)=A(T) 
END IF 
WRITE (6,600) (I,A(I),I=1,MP1) 


Calculate the Power Spectrum 


CALL SPEC (RE,IM,A,AJ,M,PI) 
D=100 
WRITE (6,601) (.005*Z,AJ(Z) ,Z=1,D) 


Plotting the Spectrum 


CALL PLOT(AJ,D) 
ISPEC=ISPEC+1000 
END IF 

CONTINUE 

WRITE (6,300) (YLL(K),K=1,NB) 
FORMAT (10(F10.7,1X)) 
FORMAT(1X,15,1X,10(F10.7,1X)) 
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600 
601 


110 


20 
200 


210 


220 


230 


240 


AAG 


260 


ae 


FORMAT (1X,15,5X,F15.7) 
FORMAT (1X,F5.3,5X,F15.7) 
STOP 
END 


SWeeOUa TNE MUMS(PH,FS,PSD,GR,GRB,GRD,F,8,BD,R,SGL,EK,AMU,N,MA) 
DIMENSION PH(100,100) ,PS(100,100) ,PSD(100,100) ,GR({100) ,GRB(100) 


1,GRD(100,100) ,F (100) ,B(100) ,BD(100) ,R(100),SGL(101) 


GR(N)=BD(N) 

GRB(N)=F(N) 

boezoo 1=2,N 

PH(I,1)=BD(N+1-I) 
PS(I,1)=F(N+1-I) 

—7= 

DO 110 K=1,1T 
PH(I,K+1)=PH(1I,K)+R(N+1-I+K)*PSD(I,K) 
PS(1I,K+1)=PSD(1,K)+R(N+1-I+K)*PH(I,K) 
CONTINUE 

DO 120 K=1,1T 

PSD(I,K)=PS(i,K) 

CONTINUE 

DO 210 K=2,N 

GR(N+1-K)=PH(K,K) 
GRB(N+1-K)=PS(K,K) 

CONTINUE 

DO 220 K=1,N 

DO 220 L=1,MA 
GRD(K,MA+2-L)=GRD(K,MA+1-L) 

DO 230 K=1,N 

GRD (K,1)=GRB(K) 

DO 240 K=1,N 

GRB(K)=GRD(K,MA+1) 

DO 250 K=1,N 

GR(K)=GR(K)+GRB(K) 

DO 260 K=1,N 
SGL(K)=.90*SGL(K)+.1LO*GR(K)*GR(K)+1.0 
BOr270 1=1,N 
R(I)=R(I)-(AMU/SGL(1)) *EK*GR(T) 
IF(R(I).GE.1.0) R(I)=0. 
IF(R(I).LE.-1.0) R(I)=0. 
CONTINUE 

RETURN 

END 


oy 


10 


Zo 


20 


a2 


oe 


oS 
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SUBROUTINE STEPUP (A,B,ML) 
DIMENSION A(100),C(100) ,B(100) 
A(1)=1. 

A(2)=B(1) 

DO 30 MINC=2,ML 

DO 10 J=1,MINC 

JB=MINC-J+1 

C(J)=A( 38) 

DO 20 IP=2,MINC 
A(IP)=A(IP)+B(MINC) *C(IP-1) 
A(MINC+1 )=B(MINC) 

CONTINUE 

RETURN 

END 


SUBROUTINE SPEC(RE,IM,A,AJ,M,PI) 

REAL RE(100),IM(100) ,A(100) ,AJ(100) , ¥(100) 
RE(1)=A(1) 

IM(1)=0. 

DO 91 J=1,100 

DO 92 I=1,M 

RE (I+1)=RE(I)+A(I+1)*COS (2*I*PI*.5*J/100) 
IM(I+1)=IM(1)+A(I+1)*SIN(2*I*PI*.5*J/100) 
CONTINUE 
AJ(J)=-10.*ALOG10(RE(M+1)*RE(M+1)+IM(M+1)*IM(M+1) ) 
CONTINUE 

TEMP=AJ(1) 

DO 93 L=2,100 

IF (AJ(L).GT.TEMP) TEMP=Ad(L) 

CONTINUE 

DO 94 L=1,100 

AJ(L)=AJ(L)-TEMP 

CONTINUE 

RETURN 

IND 


SUBROUTINE PLOT(Y,N) 
DIMENSION Y(N),X(100) 
ISTP=N/10 

DO 10 J=1,N 

X(J)=J 

DF=.5/N 
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10 


20 


X(1)= 


0 


DO 10 K=2,N 


X(K)= 


X(K-1)+DF 


XMIN=X(1) 


AMA X= 
XSTP= 


X(N) 
1O*DF 


IYMIN=Y(1L) 

IYMAX=Y (1) 

DO 20 K=2,N 
IF(¥(K).GT.IYMAX) IYMAX=Y(K) 
IF(Y¥(K).LT.IYMIN) IYMIN=Y(X) 
CONTINUE 
IYSTP=(IYMAX-IYMIN) /5 


CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 
CALL 


TEK618 

PRTPLT(72,6) 

SHERPA('PARKPARK','A' ,3) 

RESET('ALL') 

PAGE(8.5,11.0) 

HWROT('AUTO' ) 

XINTAX 

RREA2D(5.0,2.3) 

HEIGHT(0.10) 

COMPLX 

SHDCHR(90.0,1,0.002,1) 

HEADIN(' FREQUENCY SPECTRUMS! ,100,2.0,1) 
XNAME (' FREQUENCYS! ,100) 

YNAME ( 'MAGNITUDE(DB)$',100) 
MESSAG('FIGURE 5.21 (M=8 ,SNR=30DB)$',100,3.0,-0.8) 
MESSAG('MODEL ORDER SELECTION(4 SINUSOIDS)$',100,3.0,-0.8) 
THKFRM(0.03) 

FRAME 

GRAF (XMIN , XSTP, XMAX, IYMIN, IYSTP, IYMAX) 
THKCRV (0.02) 

CURVE(X,Y,N,0) 

ENDPL(0) 

DONEPL 


RETURN 


2ND 
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